###########


# Libraries
library(foreign)
library(tseries)
library(plm)
library(AER)
library(zoo)
library(TSA)
library(pcse)
library(xtable)
library(lattice)
library(Zelig)
library(lmtest)
library(nlme)
library(mvtnorm)
library(graphics)
library(Matrix)
library(apsrtable)
library(arm)

## MI Data 1

# Convert to TSCS format & attach
tsdata1 <- pdata.frame(data1, index=c("state", "year"), drop.index=FALSE, row.names=TRUE)




#### Data Management
### Computing lagged and first-differenced variables

  # Function to create lags
lag.fn <- function(x) {lag(x, 0:1)[,2]}
  # Apply function
tsdata1$L.i_entitlement <- lag.fn(tsdata1$i_entitlement)
tsdata1$L.i_socialcorp3 <-  lag.fn(tsdata1$i_socialcorp3)
tsdata1$L.leftgov <-  lag.fn(tsdata1$leftgov)

tsdata1$L.netimmrate <- lag.fn(tsdata1$netimmrate)
tsdata1$L.laborcost_new <- lag.fn(tsdata1$laborcost_new)
tsdata1$L.unemp <- lag.fn(tsdata1$unemp)
tsdata1$L.trade <- lag.fn(tsdata1$trade)
tsdata1$L.gdp <- lag.fn(tsdata1$gdp)
tsdata1$L.fdi <- lag.fn(tsdata1$fdi)


  # Function to create first-differences
diff.fn <- function(x) {diff(x, 0:1)[,2]}
  # Apply function
tsdata1$D.i_entitlement <- diff.fn(tsdata1$i_entitlement)
tsdata1$D.i_socialcorp3 <-  diff.fn(tsdata1$i_socialcorp3)
tsdata1$D.leftgov <-  diff.fn(tsdata1$leftgov)
    
tsdata1$D.netimmrate <- diff.fn(tsdata1$netimmrate)
tsdata1$D.laborcost_new <- diff.fn(tsdata1$laborcost_new)
tsdata1$D.unemp <- diff.fn(tsdata1$unemp)
tsdata1$D.trade <- diff.fn(tsdata1$trade)
tsdata1$D.gdp <- diff.fn(tsdata1$gdp)
tsdata1$D.fdi <- diff.fn(tsdata1$fdi)

# Interaction terms
tsdata1$D.netimmrateXlaborcost <- tsdata1$D.netimmrate * tsdata1$D.laborcost_new
tsdata1$D.netimmrateXsoccorp <-  tsdata1$D.netimmrate * tsdata1$D.i_socialcorp3
tsdata1$D.netimmrateXleftgov <- tsdata1$D.netimmrate * tsdata1$D.leftgov

tsdata1$L.netimmrateXlaborcost <- tsdata1$L.netimmrate * tsdata1$L.laborcost_new
tsdata1$L.netimmrateXsoccorp <-  tsdata1$L.netimmrate * tsdata1$L.i_socialcorp3
tsdata1$L.netimmrateXleftgov <- tsdata1$L.netimmrate * tsdata1$L.leftgov


#############################################################################################


# Old model interactions

tsdata1$Dimmrate_Dlabcost <- tsdata1$D.netimmrate * tsdata1$D.laborcost_new
tsdata1$Dimmrate_Lsocial <- tsdata1$D.netimmrate * tsdata1$L.i_socialcorp3



#############################################################################################


############
#  ECM
###########



ecm1 <- plm(D.i_entitlement ~ 1+
L.i_entitlement + L.netimmrate + L.i_socialcorp3 + L.leftgov +  L.laborcost_new 
+ L.unemp +  L.trade + L.gdp +  
L.netimmrateXlaborcost + L.netimmrateXsoccorp +
L.netimmrateXleftgov + 
D.netimmrate +  D.i_socialcorp3 + D.leftgov + 
D.laborcost_new + D.unemp +   D.trade + D.gdp +    
D.netimmrateXlaborcost + D.netimmrateXsoccorp +
D.netimmrateXleftgov, data=tsdata1)


ses1 <- as.matrix(sqrt(diag(vcovBK(ecm1))))

 
ecm1 <- lm(D.i_entitlement ~ 
L.i_entitlement + L.netimmrate + L.i_socialcorp3 + L.leftgov +  L.laborcost_new 
+ L.unemp +  L.trade + L.gdp +  
L.netimmrateXlaborcost + L.netimmrateXsoccorp +
L.netimmrateXleftgov + 
D.netimmrate +  D.i_socialcorp3 + D.leftgov + 
D.laborcost_new + D.unemp +   D.trade + D.gdp +   
D.netimmrateXlaborcost + D.netimmrateXsoccorp +
D.netimmrateXleftgov, data=tsdata1)

betas1 <-  as.matrix(coef(ecm1))
sigma1 <- vcov(ecm1)
est1 <- cbind(betas1[-1,], ses1)
colnames(est1) <- c("Estimates", "PCSEs")
xtable(est1, digits=3, , caption="Error Correction Model Estimates")

apsrtable(ecm1, digits=3,  coef.rows=1, caption="Error Correction Model Estimates")

#####


